package com.salinity.kun.algorithm;

public class FourierSeries {
	private double DeltaT; // 采样时间间隔
	private Double[] Ft; // 采样数据
	private int SeriesCount;// 傅里叶分级级数，默认设为T0/DeltaT/2
	private double T0; // 周期
	private int Times; // 采样点数所在时间段的总长度为周期T0的多少倍
	private double W; // 波形的角速度，值为2PI/T0

	private double[] a;
	private double[] b;

	/**
	 * 
	 * @param deltaT
	 *            //采样时间间隔
	 * @param ft
	 *            //采样数据
	 * @param seriesCount
	 *            //傅里叶分级级数，默认设为T0/DeltaT/2
	 * @param t0
	 *            //周期
	 */
	public FourierSeries(double deltaT, Double[] ft, int seriesCount, double t0) {
		DeltaT = deltaT;
		Ft = ft;
		SeriesCount = seriesCount;
		T0 = t0;
		Times = (int) (Ft.length * DeltaT / T0);
		W = 2 * Math.PI / T0;
		if (Times <= 0) {
			Times = 1;
		}
	}

	public FourierSeries(double deltaT, Double[] ft, double t0) {
		DeltaT = deltaT;
		Ft = ft;
		T0 = t0;
		SeriesCount = (int) (T0 / DeltaT / 2);
		Times = (int) (Ft.length * DeltaT / T0);
		W = 2 * Math.PI / T0;
		if (Times <= 0) {
			Times = 1;
		}
	}

	/**
	 * 获取指定index的ai值,积分单元用的是矩形，精度低于梯形
	 * 
	 * @param index
	 * @return
	 */
	public double get_a_rect(int index) {
		double a_i = 0.0;
		for (int i = 0; i < Ft.length; i++) {
			a_i += 2 * (Ft[i] * Math.cos(index * W * (i * DeltaT))) * DeltaT / (Times * T0);
		}
		return a_i;
	}

	/**
	 * [0,T]
	 * 获取 序列a，积分单元用的是梯形，精度较矩形高,公式： ai = (2/T0)∫(T0)(0) ft(i) * COS(iwt) dt
	 * 
	 * @param index
	 * @return
	 */
	public double get_a_back(int index) {
		double a_i = 0.0;
		double up, down;
		int len = Ft.length;
		
		up = Ft[0] * Math.cos(index * W * ((0) * DeltaT));
		down=Ft[len-1] * Math.cos(index * W * ((len-1) * DeltaT));
		a_i = 0.5*(up+down)*DeltaT;
		for (int i = 1; i < len-1; i++) {
			a_i+=Ft[i] * Math.cos(index * W * (i * DeltaT))*DeltaT;
		}
		return 2*a_i/(Times * T0);
	}
	
	/**
	 * [-T/2,T/2]
	 * @param index
	 * @return
	 */
	public double get_a(int index) {
		double a_i = 0.0;
		double up, down;
		int len = Ft.length;
		double T = T0;
		double T1 = T / 2f;
		double t0 = -T1;
		
		//梯形面积求和
		up = Ft[0] * Math.cos(index * W * t0);
		down = Ft[len - 1] * Math.cos(index * W * T1);
		a_i = 0.5 * (up + down) * DeltaT;
		for (int i = 1; i < len - 1; i++) {
			t0 = -T1 + i * DeltaT;
			a_i += Ft[i] * Math.cos(index * W * t0) * DeltaT;
		}
		return 2 * a_i / T;
	}

	/**
	 * 获取 序列
	 * 
	 * @return
	 */
	public double[] get_a() {
		if (a == null) {
			a = new double[SeriesCount+1];
			for (int i = 0; i < SeriesCount+1; i++) {
				a[i] = get_a(i);
			}
		}
		return a;
	}

	/**
	 * 获取指定index的 bi值，积分单元用的是矩形，精度低于梯形
	 * 
	 * @param index
	 * @return
	 */
	public double get_b_rect(int index) {
		double b_i = 0.0;
		for (int i = 0; i < Ft.length; i++) {
			b_i += 2 * (Ft[i] * Math.sin(index * W * (i * DeltaT))) * DeltaT / (Times * T0);
		}
		return b_i;
	}

	/**
	 * [0,T]
	 * 积分单元用的是梯形，精度较矩形高,公式：
	 * 
	 * bi = (2/T0)∫(T0)(0) ft(i) * SIN(iwt) dt
	 * 
	 * @param index
	 * @return
	 */
	public double get_b_back(int index) {
		double b_i = 0.0;
		double up, down;
		int len = Ft.length;
		if(index==0) return 0.0;
//		for (int i = 1; i <len ; i++) {
//			up = Ft[i - 1] * Math.sin(index * W * ((i - 1) * DeltaT));
//			down = Ft[i] * Math.sin(index * W * (i * DeltaT));
//			// 2 * 梯形面积 / T ，本来为2 * (up + down) * DeltaT / (2 * Times * T0) ，这里化简
//			b_i += (up + down) * DeltaT / (Times * T0);
//		}
		//梯形求面积另外一种方法
		up = Ft[0] * Math.sin(index * W * ((0) * DeltaT));
		down=Ft[len-1] * Math.sin(index * W * ((len-1) * DeltaT));
		b_i = 0.5*(up+down)*DeltaT;
		for (int i = 1; i < len-1; i++) {
			b_i+=Ft[i] * Math.sin(index * W * (i * DeltaT))*DeltaT;
		}
		return 2*b_i/(Times * T0);
	}


	/**
	 * [-T/2,T/2]
	 * @param index
	 * @return
	 */
	public double get_b(int index) {
		double b_i = 0.0;
		double up, down;
		int len = Ft.length;
		double T = T0;
		double T1 = T / 2f;
		double t0 = -T1;
		if(index == 0) return 0.0;
		//梯形面积求和
		up = Ft[0] * Math.sin(index * W * t0);
		down = Ft[len - 1] * Math.sin(index * W * T1);
		b_i = 0.5 * (up + down) * DeltaT;
		for (int i = 1; i < len - 1; i++) {
			t0 = -T1 + i * DeltaT;
			b_i += Ft[i] * Math.sin(index * W * t0) * DeltaT;
		}
		return 2 * b_i / T;
	}
	
	/**
	 * 获取序列b
	 * 
	 * @return
	 */
	public double[] get_b() {
		if (b == null) {
			b = new double[SeriesCount+1];
			for (int i = 0; i < SeriesCount+1; i++) {
				b[i] = get_b(i);
			}
		}
		return b;
	}

	/**
	 * 获取指定级数索引下的傅里叶展开系数值
	 * 
	 * @param index
	 * @return
	 */

	public double getFs(int index) {
		double a = get_a(index);
		double b = get_b(index);
		return Math.sqrt(a * a + b * b);
	}

	/**
	 * [0,T]
	 * 通过展开得到的系数计算F(t)，用于验证
	 * 
	 * @param index
	 * @return
	 */
	public double getFtCal_back(int index) {
		
		
		if(index <0) {
			return 0.0;
		}
		if (a == null) {
			a = get_a();
		}
		if (b == null) {
			b = get_b();
		}
		double ft = 0.0;
		ft += a[0] / 2;
		double t0 = index * DeltaT;
		
		for (int i = 1; i < SeriesCount+1; i++) {
			ft += a[i] * Math.cos(i * W * t0) + b[i] * Math.sin(i * W * t0);
		}
		return ft;
	}
	
	/**
	 * [T/2,T/2]
	 * 通过展开式计算得到deltaT整数倍时间的函数值
	 * @param index
	 * @return
	 */
	public double getFtCalByIndex(int index) {
		
		
		if(index <0) {
			return 0.0;
		}
		if (a == null) {
			a = get_a();
		}
		if (b == null) {
			b = get_b();
		}
		double ft = 0.0;
		ft += a[0] / 2f;
		
		double T = T0;
		W = 2 * Math.PI / T;
		double T1 = T / 2f;
		
		double t0 = -T1 + index * DeltaT;
		for (int i = 1; i < SeriesCount+1; i++) {		
			ft += a[i] * Math.cos(i * W * t0) + b[i] * Math.sin(i * W * t0);
		}
		return ft;
	}
	
/**
 * 
 * @param t
 *  时间t, 范围在[-T/2，T/2]
 * @return
 */
public double getFtCalByTime(double t) {
		
		if (a == null) {
			a = get_a();
		}
		if (b == null) {
			b = get_b();
		}
		double ft = 0.0;
		ft += a[0] / 2f;
		
		double T = T0;
		W = 2 * Math.PI / T;
		double t0 = t;
		for (int i = 1; i < SeriesCount+1; i++) {		
			ft += a[i] * Math.cos(i * W * t0) + b[i] * Math.sin(i * W * t0);
		}
		return ft;
	}
	
	public int getFtLength() {
		return Ft.length;
	}

}
